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We construct the S-wave part of the electromagnetic vector annihilation current to ff{asV^), where 
V is the non-relativistic quark velocity, for heavy quarks whose dynamics are described by the 
NRQCD action on the lattice. The NRQCD vector current for QQ annihilation is expressed as a 
linear combination of lattice operators with quantum numbers L = Q,J^^ 1 ^ , and the coefficients 
are determined by matching to the corresponding continuum current in QCD to ^(v^) at one- 
loop. The annihilation channel gives a complex amplitude with Coulomb-exchange and infrared 
singularities, making a careful choice for the contours of integration and infrared subtraction 
functions in the numerical integration necessary. An automated vertex generation program written 
in Python is employed, allowing us to use a realistic NRQCD action and an improved gluon lattice 
action; a change in the definition of either action is easily accommodated in this procedure. The 
final result is applicable to simulations of electromagnetic decays of heavy quarkonia, notably the 
T meson. 
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1. Introduction 

Leptonic widths of heavy quarkonia such as the T or the 7/ 1//^ are an important test of elec- 
troweak Standard Model in the heavy quark sector: The heavy quarks are the heaviest Standard 
Model particles and hence should be sensitive to possible new physics at or above the electroweak 
scale, and leptonic decays have experimentally clean signatures. Moreover, ratios of leptonic 
widths can be measured to good accuracy both experimentally and on the lattice. 

The CLEO collaboration has experimental results to few -percent precision [|l|] : 

^I^^}=flf-=QA51{6) (1.1) 

rr(15)^e+e- 



which has to be compared with the current best lattice result [g] 

-0.48(5) (1.2) 



rT(25)^e+e-^|(2S) 
rT(15)^e+e-^T(lS) 



There is thus a challenge to the lattice community to obtain a precision on theoretical predic- 
tions that can be compared to that achieved experimentally. 

2. Matching S-wave decays between NRQCD and QCD 

The leptonic width of a QQ state is given by 

^QQ^m- = KOIJ^'^'^I e!2>|'4«i (2-1) 

iMQQ 

with all the nonperturbative QCD contributions encapsulated in the matrix element (O QQ). 
Unfortunately, it is not possible to simulate QCD with heavy quarks directly due to their short 
Compton wavelengths, so Non-Relativistic QCD (NRQCD) has to be used in lattice simulations of 
heavy quarkonia. 

Hence, we need to match the desired QCD matrix element to a series of NRQCD matrix 
elements which can be measured on the lattice: 



<o|je^^|(2e> = I«,-(o 



jNRQCD 



QQ) (2.2) 



where the a, are the matching coefficients which we need to determine. For the case of S-wave 
decays, which we will study exclusively in this paper, we can take the NRQCD currents to be 

To compute the matching coefficients perturbatively, we expand both the coefficients and the 
matrix elements perturbatively: 

a,- = X (0 1 J| QQ) = Z < (0 1 J! QQf'' (2-3) 

n n 

and match order by order in a,. 
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In the T system, the order of the NRQCD expansion parameters is ^ ~ 10%. Prima 
facie, this would suggest that to achieve ^ 1% accuracy, we would need to go to ^(a^, ttsV^, v^). 
However, in matrix element ratios the ^{(xj) terms cancel, and hence we only need to include 
ff{cCs, o;vV^,v'*) corrections for ~ 1% accuracy. 

If we are only interested in the ratio of leptonic widths of, say, T{2s) and T(l5), we do not care 
about the overall normaUsation of the matrix element, and so for each decay we need only consider 
instead the quantity 

^ = (/o) + ^(/i) + ^(/2). (2.4) 
ao uq ao 

and we can define matching coefficients for the ratio as 



(0) 

" "0 "0 



(0) (1) 
CI I CZq 

.(0) 



(0) 

^2 = — = -ToY • (2.5) 

In the following, we work in the Breit frame, where the decaying meson is stationary and the quark 
has momentum = {iE, 0,0, Mv), use v as the non-relativistic expansion parameter (which is 
exact at the order to which we are working) and treat the quarks as being exactly on-shell (which 
can also be shown to be justified). 

3. The improved NRQCD action 

The improved NRQCD action used for simulations is 



with 

aHo 



A(2) 
2M 



a5H = -ci^ ^ + C2 ^(V-£-£:-V)-C3 ^a-(VxE-ExV) 

1 . aW _ (A(2))2 

2pf) ° 2M^) ~ xen{aMf 

where n is a stabiUty parameter for the euclidean-space Schrodinger equation, which must satisfy 
n>l)j {Ma) for numerical stabihty. To the perturbative order considered here, we can take c,- = 1 . 
As our glue action, we use a Symanzik improved action with tadpole improved Unks. 

4. Automatically generating Feynman rules 

In order to correctly determine the desired matching coefficients, we need to consider exactly 
the same NRQCD action in perturbation theory as is used in simulations. For the improved NRQCD 
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action, the Feynman rules become extremely complicated: The QQg vertex, e.g., has ~ 8,000 
terms, and the QQgg vertex has ~ 70,000 terms! It is clear that a traditional manual treatment 
would be extremely cumbersome and error-prone. 

For this reason, we have developed HiPPy, an automated tool for generating Feynman rules 
from lattice actions. HiPPy is written entirely in Python with companion modules in Fortran 95, 
and is freely available from any of the authors. The main strength of HiPPy lies in its great flexi- 
bility: HiPPy is capable of handling not only various kinds of NRQCD actions, but also relativistic 
(staggered, Wilson . . . ) quark and gluon actions with or without improvement. A description of 
HiPPy has been published in [^], and it is currently being used by HPQCD member to calculate 
a variety of different improvement and renormalisation constants. Due to its flexible design, a 
HiPPy-based program can easily accommodate a change in the quark or gluon action being used 
without the need for changes to the user code. 

5. Matching at tree level 

At tree-level, the relevant matrix elements are given by 



(0|jec^|(2e: 



.NRQCD 



QQ 



(0) 



(0) 



v(-p)yM(p) = x'^a + ^ Y 



where 



go{v) = 1 

g2{v) 



[May 
4 



sm 



{Maf 



aMv 



' 9 / aMv 

4sm^ 

V 2 



■ sin^(aMv) 



Expanding these matrix elements in powers of v^, we determine af^ to match: 



(0) _ 1 (0) _ 1 {aMf 



72 



(5.1) 



6. Matching to one-loop order 



Expanding the matching condition to first order in gives 
known functions of v 



(0 

wanted 



.NRQCD 



QQ 



(0) 



1 



QCD\ 



QQ) 



IQCD 



jNRQCD 



QQ 



(1) 



(6.1) 



Inrqcd 



Both the QCD and the NRQCD matrix elements on the right-hand side contain odd powers of v 
coming from the Coulomb-exchange singularity; however, only even powers of v are available for 
matching on the left-hand side, so the odd powers must cancel exactly. 
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In fact, the odd powers of v are a purely infrared phenomenon, and are known exactly: 



where h{v) is a known even function of v. We can hence analytically subtract the odd powers from 
both QCD and NRQCD by rearranging the right-hand side as 

IqCD — InRQCD = {IqCD — lodd) - i^NRQCD — lodd) \^ + hddW'i\SS (6.3) 

where ^ signifies integration over the Brillouin zone only. 

The term {Iqcd — lodd) is known analytically, while the other terms are calculated numerically 
using farmed VEGAS on the CCHPCF SunFire Galaxy class computer. 

The results obtained at various v are then fitted with the ansatz 

Iqcd - Inrqcd = a^o^ - ai^i (v) (6.4) 
to obtain the matching coefficients at one-loop order. 

6.1 The QCD form factors 

The relevant QCD on-shell form factors 

F[%AE^)f,{v^)+F^\AE^)f2{v^))x'a^f 



are simply the corresponding QED form factors (times a colour factor), which are well known at 
the one loop level. The F2 term is both UV- and IR-finite; the Fi term is UV-finite by virtue of 
the Ward identity, but IR divergences (which cancel against those in the NRQCD matrix elements) 
remain. Moreover, the F\ term contains odd powers of v which arise from the 1/v Coulomb- 
exchange singularity. As mentioned before, these odd powers are known to be the same in QCD 
and NRQCD, allowing them to be subtracted analytically. 

6.2 The NRQCD self-energy 

To account for the wavefunction renormahsation in NRQCD, as well as to establish the con- 
nection between the renormalised mass in terms of which the QCD form factors are formulated and 
the bare mass appearing in the NRQCD action, we need to compute the self-energy of the NRQCD 

heavy quarks. 

The NRQCD self-energy, which is the sum of the usual "rainbow" and "tadpole" diagrams, 
can be decomposed as 

aL{pQ,p) =A + B{po,p)aTip) + Cipo,p) [l - e"'"^" (1 - flr(p))] (6.5) 

where T{p) is the tree-level kinetic energy, and from tliis form it is straightforward to derive the 
needed quantities, namely the wavefunction and (kinetic) mass renormahsation constants as 

daT. \ 



Z^{p) = l + as(ai: + 



1 + a^lM 



d{iapo) J on-shell 
doL 



p2=0 
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Figure 1: The position of the quark (red) and gluon (blue) poles in the Minkowskified energy plane. A 
normal Wick rotation (a) is only possible for > — 2p • k; otherwise, the integration contour has to be 
deformed as in (b); for numerical work, the equivalent contour shown in (c) is adopted instead. 



Given the complicated nature of the Feynman rules, we employ automatic differentiation tech- 
niques [Q] to calculate the derivatives. 

6.3 The NRQCD vertex correction 

The NRQCD vertex correction suffers from the same infrared divergences that appear in the 
corresponding QCD diagram; we use a small gluon mass /i as an infrared regulator. 

In terms of the Minkowskified energy variable kJ^, the poles of the integrand (in the continuum 
limit) are at ±/c^ = y^k^ + jx^ - ie for the gluons, and = (^^^^M^^ ~ quarks, as 

shown in fig. |l|. Hence, a normal Wick rotation between Euclidean and Minkowskian momenta is 
possible only for > — 2p ■ k, since otherwise, the quark poles cross imaginary axis. We therefore 
need to deform the Euclidean contour of integration to avoid the quark poles and pick up the correct 
analytic continuation to Minkowksi space, and the choice of contour is shown in fig. |l](c). 

To subtract the odd powers of v coming from the Coulomb-exchange singularity, we use the in- 



tegral form of eqn. |6^. The evaluation of the resulting finite integral is still quite hard numerically, 
and takes up the major part of the computer time used. 

For / > 0, the matrix elements of the NRQCD current J, also contains tadpole-type diagrams. 
Since each tadpole loop reduces the v-dependence of the result by one power of v^, this leads to a 
contamination of the lower-order matching coefficients by "mixing down", which would appear to 
necessitate a complete recalculation if higher-order currents are added in later. The solution lies in 
defining subtracted currents Ji = ZijjJ^^'^^ where zij is defined such that we have (O |^| QQ)^"^ = 
^(v^') for all n. For details, we have to refer the reader to [S]. 



7. Results 



We have run our calculation at a number of different quark masses corresponding to the bottom 
quark on the MILC supercoarse, coarse and fine ensembles, and to the charm quark on the super- 
coarse ensemble. We have also performed extensive tests of gauge invariance, infrared regulator 
independence, and agreement with known results for at v = in the case of simpler NRQCD 
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Figure 2: Left: the numerical results with the fit to extract the matching coefficients; right: a plot of results 
in different gauges against the infrared gluon mass, showing gauge and gluon mass independence. 



Mqo n a, 



'1 



b 



4.0 2 -0.1288(27) -3.32(29) -3.30(30) -0.0972 

2.8 2 -0.1732(21) -1.35(22) -1.32(22) 0.0161 

1.95 2 -0.1358(16) 0.26(17) 0.14(17) 0.0722 

1.0 4 0.4056(20) -0.50(17) -0.56(17) 0.1111 



Table 1: The matching coefficients, as a function of the bare heavy quark mass, for the leptonic width (af) 
and leptonic w 
mixing down. 



and leptonic width ratio (Z?,). Note that Aq"' = 1, af^ = bf^ = ^, and that there is no subtraction to prevent 



actions. A plot of our results can be seen in fig. g, as can be a plot showing the gauge and regulator 
independence of our results. Our final results for the matching coefficients are given in tab. [H 
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